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Abstract. Two dimensional nuclear magnetic resonance (NMR) spectroscopy is one 
of the major tools for analysing the chemical structure of organic molecules and 
proteins. Despite its power, this technique requires long measurement times, which, 
particularly in the recently emerging diamond based single molecule NMR, limits 
its application to stable samples. Here we demonstrate a method which allows to 
obtain the spectrum by collecting only a small fraction of the experimental data. Our 
method is based on matrix completion which can recover the full spectral information 
from randomly sampled data points. We confirm experimentally the applicability of 
this technique by performing two dimensional electron spin echo envelope modulation 
(ESEEM) experiments on a two spin system consisting of a single nitrogen vacancy 
(NV) centre in diamond coupled to a single ^^C nuclear spin. We show that the main 
peaks in the spectrum can be obtained with only 10 % of the total number of the data 
points. We believe that our results reported here can find an application in all types 
of two dimensional spectroscopy, as long as the measured matrices have a low rank. 
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1. Introduction 

A key tool in the quest for the determination of the structure of molecules and 
proteins is nuclear magnetic resonance spectroscopy (NMR) which has helped to make 
fundamental contributions to the advancement of biological sciences. This is achieved by 
measuring the magnetic response of molecules in a large ensemble to sequences of radio 
frequency pulses. This temporal response is then mapped to multi-dimensional spectra 
which encode the dynamical properties of the system and therefore the interactions 
between its constituent nuclear spins [1, 2]. The information contained in these spectra 
forms the basis for the determination of molecular structure. Current NMR schemes 
are intrinsically ensemble measurements, both due to the minute size of the nuclear 
magnetic moments and the tiny polarization of these nuclear spins at room temperature, 
even in very strong magnetic helds. Consequently, NMR can only deliver ensemble 
information while the structure and dynamics of individual specimens remain hidden 
from observation. 

Recent progress in the control of the single electron spin in Nitrogen-vacancy (NV) 
centers in diamond offers a new perspective here, as it can make use of optically 
detected magnetic resonance [3, 4] for the detection of material properties [5] including 
minute magnetic helds [6, 7, 8, 9]. Building on this, recent theoretical investigations 
[10, 11, 12, 13] have suggested that NV centers implanted a few nanometers below the 
surface should be able to detect and locate individual nuclear spins above the diamond 
surface. Subsequent experimental work has indeed achieved the observation of small 
clusters of nuclear spins outside of diamond with a sensitivity that is sufficient to identify 
even individual nuclear spins [14]. 

One of the challenges for the determination of the structure of smaller biomolecules 
or even entire proteins by means of 2D spectroscopy detected by an NV center is the 
considerable amount of data that need to be taken which results in long measurement 
times. Indeed, the large amount of required data and the associated long measurement 
times represent a challenge that is common to both ensemble NMR and single molecule 
NMR measurements. 

As suggested in [12] we demonstrate NV sensing experiments on nuclear spins using 
methods from the held of signal processing, particularly matrix completion [15, 16]. 
With this technique we can obtain reliably the spectral information that is contained in 
2D-NMR spectra from a small subset of all accessible data points (see [17] and [18] for 
applications of the related but distinct compressive sensing and non-uniform sampling 
to bulk NMR). The results presented here show that order of magnitude reduction in 
the overall measurement time in NV center based 2D-NMR can be achieved. 

In the remainder we briehy introduce matrix completion in Section 2. Then Section 
3 presents the application of this method to concrete experimental data that have 
been obtained from an NV center interacting with a nearby nucleus. The results 
demonstrate that already a sampling rate of around 10% suffices to reconstruct the 
spectral information reliably. We hnish with a brief conclusion and outlook concerning 
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2. Matrix Completion Method 

This section serves to introduce briefly the concept of matrix completion, the basic 
properties relevant to this work and the specihc algorithm that we use for its application 
to our experimental data. 

A 2D-spectrum encodes the response of a system to a sequence of pulses with 
varying temporal separation, denoted by ti and Tj, and arrange the result in a data 
matrix M{i,j). The 2D-spectrum S'(a;i,a; 2 ) is then obtained as the Fourier transform 
of both time coordinates in M. In our work we are sampling randomly chosen elements 
of the matrix M with indices (f,j) drawn from the index set leading to constraints 
Xij = Mij for {i,j) G Matrix completion solves the task of obtaining the missing 
matrix entries of M that have not been measured in experiment. In general this is 
impossible unless we have further knowledge about the matrix M, namely that it 
typically has a low singular value rank r, i.e. r n. Fortunately, this is indeed 

the case for typical data sets and especially those from 2D spectroscopy. 

One possible approach to achieve this matrix completion is by solving the 
minimization problem 

min [tr|X| : \Xij — Mij\ < e for (z,j) G D] (1) 

where tr|X| is the trace norm of the matrix X and e is a given tolerance. Indeed, it can be 
proven that this formulation of the problem achieves the desired aim [19] as the solution 
of eq. (1) yields the matrix M with high probability if the number of sampled elements 
|f2| = 0{nr\\in) where r is the singular value rank of the n x n matrix M (see [19, 20] 
for proofs and a rigorous mathematical statement). This suggests that a computational 
gain by a factor of order n/{rlnn) may be achieved through random sampling in the 
manner described above (see for example [12, 21] on computed 2D-spectroscopy data). 

This still leaves us with the task of solving the minimization problem eq. (1). 
In principle, this equation can be rewritten as a semi-dehnite programme and then 
solved employing standard solvers for convex problems. Unfortunately, standard solvers 
tend to be limited to relatively small matrix sizes, but fortunately alternatives exist. 
Indeed, [22] proposed to solve eq. (1) approximately through the so-called singular value 
thresholding (SVT) algorithm [22] which permits very large matrices to be treated. It is 
this algorithm that we will be using in our work. The SVT-algorithm solves iteratively 
the set of equations 

y(fc-l) _ jj{k-l)y{k-l) 1 ^ 2 ^ 

Xik) = f/(fc-i) max - rl) (3) 

Y(k) = yU-i) + 5kVn{.M - X^^'>) (4) 

where {Vn{M))ij = Mij for (hi) G D and zero otherwise and eq. (2) represents the 
singular value decomposition of the matrix ^ and 5^ are free parameters in the 
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procedure that regulate the soft thresholding (eq. 4) and the inclusion of the constraints 
(eq. 5). The choice 5^ < 2 ensure provable convergence and t = 5n for n x n-matrices 
represent typical values (see [22]). As a termination criterion of the iteration we employ 
the condition 

|(P„(X(‘I-M)|| 

ll(Pn(M)||^ ^ ' ' ' 

for some e, and ||.||f being the Frobenius norm [22], The algorithm employs a 
singular value decomposition which, for large matrices, can be accelerated considerably 
[23, 24], It is also noteworthy that other approaches for solving eq. (1) such as those 
reported in [25, 26, 27] may lead to improved performance and/or stability but for the 
purposes of this study SVT was sufficient and recommended itself thanks to its ease 
of implementation. In any real-world application, the measured entries of the data 
matrix will be corrupted at least by a small amount of noise. Hence the question of 
the robustness of the matrix completion approach to fluctuations in the experimental 
data arises naturally. Reassuringly, results have been developed that guarantee that 
reasonably accurate matrix completion is possible from noisy sampled entries [28]. 
In that scenario noise can be neglected if the relevant spectral information can be 
still extracted from the low rank approximation of M, thus implicating a sufficiently 
large signal-to-noise ratio and leading to the fact that noise contribution results in 
small singular values, which are discarded after applying our algorithm. Hence matrix 
completion offers three major advantages: 

• Weak noise is directly suppressed by the matrix completion algorithm 

• The spectrum of the system can be recovered from a small subset of all data e.g. 
only 10 % of the total in our examples. 

• In contrast to compressed sensing, the algorithm used here does not require any 
additional information, e. g. the sparse basis [29]. 

The following section will now present the result of the application of the matrix 
completion algorithm to concrete experimental data that have been obtained in our 
laboratory. 

3. Experimental Implementation 

3.1. 2D ESEEM with a single NV centre 

The method of matrix completion has been implemented in 2D optical spectroscopy 
of Rb vapour [30]. We use a single NV centre in diamond coupled to a proximal 
nuclear spin as a test system for the demonstration of the matrix completion protocol. 
NVs are optically active point defect centres in the diamond crystal. Their fluorescence 
depends on the electron spin number of the triplet ground state, allowing to measure 
the electron spin of single centres. NVs close to the diamond surface have been used 
to detect few thousand external protons [31, 32] followed later by a demonstration 
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of even single spins sensitivity [14] leading to nano-scale magnetic resonance imaging 
[33, 34, 35]. For these type of experiments the data acquisition is quite long due to the 
low fluorescence emission from single centres. 

The NV has a triplet ground state (electron spin S' = 1) coupled to the nitrogen nuclear 
spin / = 1). The system can be described by the Hamiltonian: 

H = DS^^ + ^B-S + S-Ai4^-I ( 6 ) 

n 

where D/2tt = 2.87 GHz is the zero held splitting of the ground state, g = 2.003 is the 
Lande factor, is the Bohr magneton, B = B^Cx -l- ByCy -|- B^Cz is the applied static 
magnetic held, S = Sx + Sy + Sz and I = Ix + ly + h are the electron and nuclear spin 
operators and is the hyperhne interaction tensor. The z axis is taken to be along 
the NV crystal axis. If there is a single nuclear spin (/ = 1/2) in the proximity, the 
following term 

Hbf^sc = S ■ Ai3c ■ / (7) 

is added to the spin Hamiltonian (6), with Ai 3 c being the hyperhne interaction tensor 
to a nuclear spin. One of the simplest 2D NMR experiments consists of three 
7r/2 pulses and is called correlation spectroscopy (COSY) [36]. In our work we use its 
’’equivalent” in the electron spin resonance - the three pulse electron spin echo envelope 
modulation (ESEEM) pulse sequence (see [37] for more details) shown in hgure 1. 


Polarize Readout 



time 


Figure 1. Pulse sequence for the two dimensional ESEEM measurement used in our 
experiments. 


The sequence starts with a laser pulse of about 3 /ts to polarize the NV electron 
spin in the \ms = 0) state. Afterwards we apply four tt/2 microwave pulses at times 
f = 0, f = Ti, f = Ti -|- r 2 and t = 2 ti + T 2 . The last pulse is used to transfer the electron 
spin coherence into population, which is read out by the last laser pulse. The spin 
signal is recorded for each pair of (ri,r 2 ) and then a 2D Fourier transform is performed 
giving a set of frequencies (i/i,i/ 2 )- From this spectrum the number of nuclei coupled to 
the electron spin and the oh diagonal elements of the hyperhne interaction tensor (e.g. 
proportional to Szlx) can be obtained [37]. 

We applied this pulse sequence in two diherent experiments. In the hrst measurement 
we use a single NV without resolvable coupling to spins. The system consists of a 
NV electron and a nitrogen nuclear spins, which are described by the Hamiltonian in 
(6). If the static magnetic held is aligned with the NV axis, the hyperhne interaction 
tensor Ai4n is diagonal and there is no ESEEM ehect. In order to introduce artihcial 
”oh-diagonal” terms, we apply the static |i?| ~ 100 G oh-axis, at an angle of about 34 
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degrees with respect to the 2 ; axis. The expected 2D spectrum S'theo can be simulated by 
using the Hamiltonian (6) and is plotted in hgure 2a. In hgure 2b we plot the Fourier 


Simulation a) 






Figure 2. 2D ESEEM simulation and experimental data with a single NV when 
static magnetic field Bq = 100.9 G is applied at an angle of 34.1°. (a) Simulation using 
the Hamiltonian (6). (b) Fourier transform of the complete set of the experimental 
data points Fourier transform of the experimental data after applying matrix 

completion and using 20 % (c) and 5'exp^ 10 % (d) of the time domain data. The 

main peaks are still observed even when 90 % of the data is removed! 

transform of the experimental data where we nse all collected data points. The 
experiment agrees well with the simnlation. In order to demonstrate the performance of 
the matrix completion method, we use a random mask D to remove a certain part from 
the fnll experimental data in the time domain. After that, we apply matrix completion 
nsing the SVT algorithm as described in section 2 to recover the fnll matrix. A Fonrier 
transform of the matrix obtained with 20 % of the initial data points 3^^^° is shown in 
hgure 2c. Despite the removal of 80 % of the recorded data, the number of peaks and 
their positions are still present if we compare to hgure 2b. Even if we keep only 10 % of 
the original matrix (see hgure 2d), we can still recover the relevant spectral information. 
In the second experiment we localized an NV coupled to a single spin with a conpling 
strength of Ahq = 9 MHz. Now, depending on the position of this carbon atom with 
respect to the NV, there are diherent hyperhne interaction tensors [38, 39, 40]. The 
spectrnm can be calculated by using the Hamiltonian (6) and (7) by choosing the 
correct hyperhne interaction tensor. The simulation is shown in hgure 3a. In these 
measnrements the magnetic held has been aligned along the NV axis. The 2D Fonrier 
transform of the fnll data set is shown in hgure 3b. As in the previous experiment, we 
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Simulation a) 






Figure 3. 2D ESEEM simulation and experimental data with a single NV coupled 
to a single nuclear spin, (a) Simulation using the spin Hamiltonian (6) and (7). 
(b) Fourier transform of the complete set of the experimental data points. Fourier 
transform of 20 % (c) and 10 % (d) of the data in the time domain. Here we again can 
recover the spectral information by keeping small amount of the experimental data. 


can still recover the full spectral information (cf. figure 3c), if we remove randomly 80 % 
of the data points. From hgure 3d we can conclude that even 10 % of the data suffice 
for the matrix completion algorithm to obtain the spectrum. In fact, this factor of ten 
is what is expected from the theory, see section 2 and below. 


3.2. Performance of the matrix completion algorithm 


In the following the performance of the matrix completion algorithm will be analysed 
on our experimental data. For this purpose we have to quantify how good we can 
recover the matrix containing the spectrum, when a small number of measurements is 
performed. We use the experimental data shown in hgure 3. The data is stored in a 
matrix Mred with reduced number of elements different from zero. This matrix is then 
compared with the matrix Mtot, obtained by measuring the total number of points with 
all pairs of (rj,rj). By using these two matrices, we dehne the hdelity F of our algorithm 
as 


F = 1- 


||M, 


tot 


M 


red 11 JP 


( 8 ) 

In hgure 4 we plot the hdelity as a function of the fraction of the elements of the 
complete matrix (red markers). In the same plot we show the ordered singular values 
(Tj of the matrix with the full number of points Mtot- From the plot we can conclude 
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Figure 4. Mean of fidelities (red dot markers, right axis) of the matrix completion 
algorithm (r = 100) as a function of the fraction of sampled elements (top axis) 
taken from the experimental data shown in figure 3. The matrix completion algorithm 
was performed 128 times with each time different random sampling, here the errorbars 
denote the standard deviation. The blue curve (left axis) represents the singular values 
CTi of the measured data (fig. 3b) in descending order. 


that most of the spectral information (> 70%) can be recovered by using 10% of the 
elements of Mtot, since only these elements are signihcantly larger than zero. This result 
is consistent with the theoretical limit for recovering Mtot given by nr Inn [19, 20] where 
we can roughly assume r = 4, which is the number of the peaks in the spectrum. The 
spectrum consists of few peaks, while the rest of the matrix elements contain only noise 
and can be discarded. 

It is interesting to investigate the influence of the threshold parameter r on the 
performance of the SVT algorithm and the hdelity of the so determined spectra (see 
Figure 5). Too small threshold values r, e.g. at r = 10 or r = 20 (pink and orange 
markers), lead to low hdelity when less than 60 % of the matrix elements are sampled. 
We can achieve higher hdelities F by increasing r and we observe saturation around 
r = 100. That is, for r ^ 100 we cannot obtain higher signihcantly higher hdelities, 
while the required computation time (equivalent to the number of iterations) increases 
which can be seen in the inset graph in hgure 5. From there we can conclude that with 
our data set thresholds even below the empirically suggested rule r = bn k, 1000 for 
our case of n = 201 (see section 2 and [22]) are sufficient. A python script can be found 
in the supplementary information, where the SVT algorithm is implemented, together 
with the data set from hgure 3b. 
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Figure 5. The fidelity of the matrix completion algorithm as a function of the fraction 
of the matrix elements at different thresholds r. Inset: Number of iterations required 
to run the matrix completion algorithm as a function of the threshold and the fraction 
of the matrix elements. The performance calculation are based on data shown in 
figure 3. 


4. Conclusions 

In summary, we have demonstrated the application of a method for reconstructing a 
two dimensional ESEEM spectrum, by collecting only small part of the data in the 
time domain. With our technique we can obtain the necessary spectral information by 
measuring 10 % of the experimental data points in two different experiments. By using 
our method, the measurement time can be shortened by a factor 10 at the same signal- 
to-noise ratio compared to the conventional experiment. We believe that the reported 
results will be useful for any type 2D NMR and ESR spectroscopy and also for magnetic 
resonance imaging. Our method is particularly useful for single spins experiments, which 
usually require very long measurement times [35, 13]. 
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